home *** CD-ROM | disk | FTP | other *** search
/ Shareware Grab Bag / Shareware Grab Bag.iso / 090 / byt1186b.arc / APPROX.LBR / CHEBY.C < prev    next >
Text File  |  1986-04-11  |  3KB  |  150 lines
  1. /*    cheby.c
  2.  *
  3.  * Program to calculate coefficients of the Chebyshev polynomial
  4.  * expansion of a given input function.  The algorithm computes
  5.  * the discrete Fourier cosine transform of the function evaluated
  6.  * at unevenly spaced points.  Library routine chbevl.c uses the
  7.  * coefficients to calculate an approximate value of the original
  8.  * function.
  9.  *    -- S. L. Moshier
  10.  */
  11.  
  12. extern double PI;        /* 3.14159...    */
  13. extern double PIO2;
  14. double cosi[33] = {0.0,};    /* cosine array for Fourier transform    */
  15. double func[65] = {0.0,};    /* values of the function        */
  16.  
  17. main()
  18. {
  19. double c, r, s, t, x, y, z, temp;
  20. double low, high, dtemp;
  21. long n;
  22. int i, ii, j, n2, k, rr, invflg;
  23. short *p;
  24. char st[40];
  25. double cos(), log(), exp(), sqrt();
  26.  
  27. low = 0.0;        /* low end of approximation interval        */
  28. high = 1.0;        /* high end                    */
  29. invflg = 0;        /* set to 1 if inverted interval, else zero    */
  30. /* Note: inverted interval goes from 1/high to 1/low    */
  31. z = 0.0;
  32. n = 64;            /* will find 64 coefficients            */
  33.             /* but use only those greater than roundoff error */
  34. n2 = n/2;
  35. t = n;
  36. t = PI/t;
  37.  
  38. /* calculate array of cosines */
  39. puts("calculating cosines");
  40. s = 1.0;
  41. cosi[0] = 1.0;
  42. i = 1;
  43. while( i < 32 )
  44.     {
  45.     y = cos( s * t );
  46.     cosi[i] = y;
  47.     s += 1.0;
  48.     ++i;
  49.     }
  50. cosi[32] = 0.0;
  51.  
  52. /*                            cheby.c 2 */
  53.  
  54. /* calculate function at special values of the argument */
  55. puts("calculating function values");
  56. x = low;
  57. y = high;
  58. if( invflg && (low != 0.0) )
  59.     {    /* inverted interval */
  60.     temp = 1.0/x;
  61.     x = 1.0/y;
  62.     y = temp;
  63.     }
  64. r = (x + y)/2.0;
  65. printf( "center %.15E  ", r);
  66. s = (y - x)/2.0;
  67. printf( "width %.15E\n", s);
  68. i = 0;
  69. while( i < 65 )
  70.     {
  71.     if( i < n2 )
  72.         c = cosi[i];
  73.     else
  74.         c = -cosi[64-i];
  75.     temp = r + s * c;
  76. /* if inverted interval, compute function(1/x)    */
  77.     if( invflg && (temp != 0.0) )
  78.         temp = 1.0/temp;
  79.  
  80.     printf( "%.15E  ", temp );
  81.  
  82. /* insert call to function routine here:    */
  83. /**********************************/
  84.  
  85.     if( temp == 0.0 )
  86.         y = 1.0;
  87.     else
  88.         y = exp( temp * log(2.0) );
  89.  
  90. /**********************************/
  91.     func[i] = y;
  92.     printf( "%.15E\n", y );
  93.     ++i;
  94.     }
  95.  
  96. /*                            cheby.c 3 */
  97.  
  98. puts( "calculating Chebyshev coefficients");
  99. rr = 0;
  100. while( rr < 65 )
  101.     {
  102.     z = func[0]/2.0;
  103.     j = 1;
  104.     while( j < 65 )
  105.         {
  106.         k = (rr * j)/n2;
  107.         i = rr * j - n2 * k;
  108.         k &= 3;
  109.         if( k == 0 )
  110.             c = cosi[i];
  111.         if( k == 1 )
  112.             {
  113.             i = 32-i;
  114.             c = -cosi[i];
  115.             if( i == 32 )
  116.                 c = -c;
  117.             }
  118.         if( k == 2 )
  119.             {
  120.             c = -cosi[i];
  121.             }
  122.         if( k == 3 )
  123.             {
  124.             i = 32-i;
  125.             c = cosi[i];
  126.             }
  127.         if( i != 32)
  128.             {
  129.             temp = func[j];
  130.             temp = c * temp;
  131.             z += temp;
  132.             }
  133.         ++j;
  134.         }
  135.  
  136.     if( i != 32 )
  137.         {
  138.         temp /= 2.0;
  139.         z = z - temp;
  140.         }
  141.     z *= 2.0;
  142.     temp = n;
  143.     z /= temp;
  144.     dtemp = z;
  145.     ++rr;
  146.     sprintf( st, "/* %.16E */", dtemp );
  147.     puts( st );
  148.     }
  149. }
  150.